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ABSTRACT 


This report presents a comprehensive summary of a number of methods used in 
filtering narrowband jammers and interference out of direct sequence spread spectrum 
communication systems. The basic concept is presented first followed by a discussion of 
several analog techniques for interference suppression. Most of the emphasis is on digital 
adaptive methods, which rely on adaptive filtering algorithms. Many of these algorithms 
are highlighted. This is followed by a discussion on the research which has been conducted 
in digital excision over the past several years. Two cases based on the author’s work are 
also presented. The first details the performance one can expect from the recursive least 
squares algorithm when the interferer is a stable tone. A.n analytical model, which predicts 
the bit error rate performance for the recursive algorithm under certain conditions, is 
derived. This model is compared to results obtained from computer simulation. The 
second case compares, through computer simulation, the tracking and filtering capabilities 
of the block and recursive least square algorithms when the interference is of the FM—type, 
i.e., a swept tone having a bandwidth which is approximately 20% of the chip rate. 



RESUME 

Ce rapport presente un sommaire complet de plusieurs metodes de filtrage de 
brouillage et d’intererences a etalement de spectre utilisani la modulation directe par 
sequences pseudo—aleatoires. On presente d’abord le concept de base et ensuite, une 
analyse de plusieurs techniques analogues de suppression d’interference. L’accent est mis 
sur des methodes adaptable numeriques reposant sur des algorithmes de filtrage adaptable 
dont plusieurs sont decrits brievement. On poursuit par une revue de la recherche effctuee 
en excision numerique au cours des dernievement. On poursuit par une revue de la 
recherche effectuee en excision numerique au cours des dernieres annees. Deux cas tires des 
travaux de l’auteur sont presentes. Le premier cas decrit les performances que l’on peut 
esperer d’un algorithme recursif a variance minimale lorsque l’interference est une 
frequence continue stable. On en derive un modele analytique qui, sous certaines 
conditions, predit les performances en terme du taux d’erreur des bits. Ce modele est 
compare a des resultats obtenus d’une simulation par ordinateur. Pour le second cas, on 
compare a l’aide d’une simulation par ordinateur, les capacites de filtrage et de poursuite 
d’un algorithme recursif a variance minimale et d’un algorithme de bloc pour des 
interferences de type FM, i.e. une modulation de frequence en dent de scie dont la bande 
passante est d’environ 20% la frequence des sauts de phase. 
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EXECUTIVE SUMMARY 


This report presents a comprehensive summary of the research which has been 
conducted in the area of filtering narrowband interferes out of direct sequence spread 
spectrum signals. These signals are used extensively in military communication systems, 
especially satellite systems. The techniques described herein apply equally to both 
Electronic Support Measure (ESM) systems and direct sequence spread spectrum 
communication systems. In the former application, the ESM system may be attempting to 
intercept the spread spectrum signal, but the narrowband interference may be hampering 
this effort. In the latter application, the spread spectrum communication system may 
require additional processing gain to aid in suppressing the interference. Since the open 
literature has been devoted to this latter case, the material presented here focuses on this 
application. 

One of the attributes of direct sequence spread spectrum communication systems is 
their ability to combat interference or intentional jamming by virtue of the system’s 
processing gain inherent in the spreading and despreading process. The interference can be 
attenuated by a factor up to this processing gain, which may be in the range of 20 to 40 
dB. In some cases this gain is insufficient to effectively suppress the interferes leading to a 
significant degradation in communications manifested by a sudden increase in bit error 
rate. If the ratio of interference bandwidth to spread spectrum bandwidth is small, the 
interference can be filtered out to enhance system performance. However, this is at the 
expense of introducing some distortion onto the signal. This process of filtering is 
sometimes referred to as interference excision. 

There are a host of analog and digital techniques which can be used to effect 
excision. Three of the more well-known analog methods are described in this report. They 
include the use of a notch filter (which may be tunable), phase—locked loops, and surface 
acoustic wave devices. The problem typically associated with analog methods is that they 
are unable to cope with a changing interference environment. Digital techniques are more 
amenable to this situation. 

Digital techniques require that the received composite signal, consisting of the 
spread spectrum signal, noise and interference, first be digitized. The resulting signal is 
then applied to a digital filter which carries out the process of excision. Its filter 
coefficients are obtained from an algorithm usually based on a minimum mean squared 
error optimization criterion. These algorithms can be of the block or recursive types. The 
block algorithms calculate a set of filter parameters from a block of data. This same block 
of data is then applied to the filter. The recursive algorithms, on the other hand, update 
the filter coefficients every data sample. This report provides a comprehensive summary of 
both groups of algorithms. 

These algorithms have formed the basis for an extensive amount of research over the 
past several years in interference excision. The salient aspects of this research are 
prestuicd in chronological order. The interference conditions included a single tone, 
multiple tones, coloured Gaussian noise, a pulsed tone, and in one instance a swept tone. 
These include common interference and jamming signals 






Several filter configurations also evolved over time, commencing with the simple 
linear prediction and interpolation filters and branching out into non—linear filters, which 
use decision—feedback schemes. The non—linear methods provide a substantial 
improvement in performance over their linear counterparts, and approach ideal 
performance. 

The report concludes with work carried out by the author. Two cases are presented. 
The first details the performance that can be expected from the recursive least squares 
algorithm when the interferer is a stable tone. An analytical model, which predicts the 
bit-error—rate performance for the recursive algorithm under certain conditions, is derived. 
This model is compared to results obtained from computer simulation. The second case 
compares, through computer simulation, the tracking and filtering capabilities cf the block 
and recursive least squares algorithms when the interference is of the FM—type, i.e., a 
swept tone having a bandwidth whica is approximately 20% of the chip rate. 

The results for this second case show that one does not have to resort to high order 
filters to suppress the tone, a second order filter will be equally effective. This reduces the 
number of mathematical operations per iteration required for calculating the filter 
coefficients. These results contrast those found for the stable tone interferer, in which 
increasing the filter order improved the performance of the excision filter. 

Finally, for an FM-type of interferer, the recursive algorithms produce a notch 
which lags behind the interference. The amount of lag is a function of the parameters used 
in the recursive algorithm and the rate of change of the modulating signal in the FM 
interferer. These findings have led the author to suggest another scheme to filter such an 
interferer. The report concludes with the proposed concept. 

Military applications of this work cannot be discussed in this unclassified report. 
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1.0 INTRODUCTION 


Direct sequence spread spectrum (DS/SS) communication systems have an inherent 
processing gain which has the effect of reducing the harmful effects of interference. When 
the interference has a power advantage over the spread spectrum system, a severe 
degradation in communication can result as manifested by a sudden increase in the bit 
error rate (BER). The processing gain can be enhanced somewhat by filtering the 
interference, particularly if it is narrowband. This technique is sometimes referred to as 
interference excision. 

Several methods of interference excision are presented in this report. First, the 
basic concept is described. This is followed by a summary of some of the analog techniques 
to implement excision filters, which may involve the use of simple notch filters, 
phase—locked loops (PLL’s), or surface acoustic wave (SAW) devices. 

Most of the emphasis in this report is on adaptive digital filtering methods based on 
time—domain modeling techniques. There are a number of algorithms associated with these 
models, and many of these will be discussed. Following this, the wc*k of numerous 
researchers in the area of digital excision will b" summarized. The report concludes with 
work conducted by the author on the application of two adaptive algorithms (block least 
squares and recursive least squares) to the filtering of stable and swept tones from a binary 
phase—shift keyed (BPSK) spread spectrum system. 

2.0 BASIC CONCEPT 

Consider the received spread spectrum signal of bandwidth B ss = 2R C as shown in 
Fig. 1, centered at an intermediate frequency (IF) of f 0 Hz, where R c is the chip rate of the 
spread spectrum signal. Also present in the signal bandwidth is an interferer of bandwidth 
Bj << B ss , and centered at f 0 + Af. The interference can be filtered out eitlu - 
at the IF or at baseband. Two possible excision filter configurations are shown in Fig. 2. 
Fig. 2(a) illustrates an analog excisor and Fig. 2(b) illustrates a digital one. In these 
figures r(t) is defined as 

r(t) = s(t) + n(t) + i(t), (1) 

where s(t) is the spread spectrum signal, n(t) is bandpass Gaussian noise and i(t) is the 
narrowband interring signal. These three components are further defined as 


s(t) = d(t)p(t)cos(i t , 0 t) 

(2a) 

n(t) = v/ n c (t) J + n s (tpcos[i- 0 t + d(t)] 

(2b) 

i(t) = A(t)cos[t' 0 t 4- 0(t)] 

(2c) 


where d(t) is a random sequence of bits, e°ch of duration Tb seconds and amplitude (± 1), 
p(t) is a pseudo random sequence of chips, each of duration T c seconds and amplitude 
(± 1), u 0 = 2rf 0 where f 0 is the carrier, n c (t) and n s (t) are in-phase and quadrature 
baseband Gaussian noise components of variance NoB ss where N 0 is the single-sided power 
spectral density of the noise, Mi) = arctan[n s (t)/n, : (t)j is the phase angle between the 
Gaussian noise components, A(t) is the time—varying amplitude of the interferer, and 0(t) 
is the phase of the interferer. 
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FIGURE 2 (a): ANALOG EXCISOR. 


Figure 2.2 (a) 



FIGURE 2 (b): DIGITAL EXCISOR. 










In Fig. 2(a), r(t) is applied to an excision filter first. The output is despread by 
correlating it with the pseudo—noise (PN) sequence, mixed to baseband and bit-detected. 
In Fig. 2(b), r(t) is mixed to baseband first and then applied to an integrate—and—dump 
circuit whose output is sampled at the chip rate. The assumption is made that carrier and 
chip synchronization have been achieved, which will not be the case initially. The sampled 
sequence, x n , is processed by an algorithm which determines a set of filter coefficients for 
the excision filter using some optimization criterion. The output of the filter, f n , is then 
correlated with the PN sequence and applied to the bit detector. 


3.0 ANALOG EXCISION TECHNIQUES 

Several analog techniques for suppressing interference in spread spectrum systems 
have been proposed [79],[3],[52],[50]. Three of the more well-known methods will now be 
discussed. 

The simplest method [79] is to apply the received signal to a tunable notch filter of 
bandwidth sufficiently large to attenuate the interference. The n^tput of the filter will 
consist of residual interference, noise and a somewhat distorted signal due to the filtering 
operation. This latter effect degrades the correlation of the signal with the PN sequence. 
The amount of degradation depends on the desired attenuation and the filter order and 
type. 


Experimental results of these tradeoffs were presented in [79]. The interference was 
a double—sideband suppressed carrier signal occupying a bandwiath equal to 3% of the chip 
rate and centered at the carrier of the spread spectrum signal. A second order Butterworth 
filter with a 10 dB bandwidth 2.5 times larger than the interfered 10 dB bandwidth 
achieved at least 25 dB interference suppression. The penalty in this case was a 1 dB drop 
in the PN correlation peak at the output of the PN correlator shown in Fig. 2(a). For the 
same interference conditions and a fourth order Butterworth filter, the degradation in the 
correlation peak was 1.75 dB. Results also showed that the amount of degradation lessened 
as the interference and notch filter were offset from the carrier. Finally, theoretical curves 
indicated that Tchebycheff filters introduced slightly more degradation in the correlation 
peak than the Butterworth filters of the same order. 

A more sophisticated analog technique [3], illustrated in Fig. 3, uses a PLL. In this 
application, the received signal is applied to a PLL of bandwidth B^. For the condition of 

large (> 10 dB) interference—power-to-noise-plus—signal-power ratios, the PLL locks on 
to the interference. The voltage controlled oscillator (VCO) signal, which is in quadrature 
to the interference, is phase-shifted by x/2 radians and applied to two mixers Ml and M2. 
Ml mixes the received signal r(t) to baseband. The baseband signal, yft), is applied to a 
low pass filter of bandwidth B a . The output is an estimate of the interference amplitude, 

A(t), contaminated by noise and spread spectrum signal. The second mixer M2 mixes the 

same t/ 2 phase—shifted VCO signal from the PLL with A(t) to form an estimate of the 
interference, i.e., 

i (t) = A(t)cos[i' 0 t + 0(t)] (3) 
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FIGURE 3: AN ANALOG EXCISOR IMPLEMENTED WITH A PHASE-LOCKED 
LOOP. 
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which is subtracted from r(t). An analysis carried out in [3] determined that the residual 
interference—to—signal ratio at the input to the spread spectrum processor is (B^+BJ/B SS . 

This result suggests that, in theory, a significant amount of suppression of the interferer 
can be achieved, especially for stable tones in which B„ and B a can be quite small in 

comparison with the spread spectrum bandwidth, B ss . 

The third analog approach to be discussed implements excision with SAW devices. 
This method is referred to as transform domain filtering. A significant amount of research 
has been conducted in this area [9],[14],[52]—[55]. The concept is illustrated in Figs. 4(a) 
and lb). The received signal is first windowed to Tb seconds by a window function w(t) 
which, for example, can be rectangular or Hamming. The windowed signal is Fourier 
transformed by the first SAW device. The transform evolves over time, becoming valid 
when the windowed signal is fully contained in the device. The transform, therefore, is 
valid over the time interval T b to Ti, where T t is the interaction time of the SAW device. 
The output of the first device is multiplied by a transfer function, K c (4/), which is related 
to a notch filter with transfer function Hr(w), i.e., 

H c (w) = 4H R (4/)cos(i/Ti) (4) 

where the cosine term adds a time delay of T| seconds to the input signal (since the inverse 
transform calculated by the second SAW device starts to become valid at Ti seconds). The 
frequency variable u is a function of time and is equal to 2/?t, where /?, which is defined 
below, is a parameter of the SAW device. As the spectrum evolves over time, the portion 
containing the interference is "blanked out" by a gating function represented by H c (w) in 
Eq. (4). The shape of H c (t’) is rectangular. The notched spectrum is inverse transformed 
by the second SAW device. This transform produces a time domain signal of duration T b 
seconds valid over the range [Ti,Ti+T b ]. This time domain signal consists of a sequence of 
chips which are then correlated with the appropriate section of the PN sequence, shown in 
Fig. 4(a) as the spread spectrum processor. The output of this processor is then 
bit-detected. 

A detailed representation of the above process is illustrated in Fig. 4(b). Within the 
dashed lines is the way in which the Fourier transforms are implemented with the SAW 
devices. The windowed waveform is mixed with a chirp signal defined as cos(&/ a t — /?t 2 ) 
where w a is the carrier frequency in radians per second and /? is the phase sweep rate in 
radians/second 2 of the chirp. The product is applied to a tapped delay line, constructed on 
a SAW device, with an impulse response cos(t» a t 4- /?t 2 ). The inverse Fourier transform is 
similarly constructed. The PN correlation process at the output of the second SAW device 
is represented in Fig. 4(b) as a filter with impulse response s(T b — t). This is the matched 
filter (code sequence) for the particular bit, and is an alternate representation for the 
correlation process which was shown in Fig. 2(a) [83]. 

Milstein and Das [53] have used the model in Fig. 4(b) to develop a test statistic 
[83] at the output of the bit detector. This test statistic is defined as 

T b 

U = J r(t)[h R (t)*s(t)]dt (5) 

0 


- 6 - 








FIGURE 4: AN ANALOG EXCISOR IMPLEMENTED WITH SAW DEVICES; 

(a) CONCEPTUAL BLOCK DIAGRAM 

(b) IMPLEMENTATION. 
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where hR(t) = J- 1 {Hr(w)}, 7-» refers to the inverse Fourier transform and * refers to the 
convolution operator. From the mean and variance of U for a given interference and 
filtering condition, a theoretical bit error rate expression can be determined. Milstein et al. 
presented results for two such cases, i.e., a tone interferer and Gaussian noise jammer. 

They also presented an adaptive scheme with supportive experimental results [78]. 

As a perspective on the performance of such a SAW excision system, and the 
parameter values which are achievable, one can refer to the experimental system in [67]. 
The SAW devices used therein had center frequencies of t/ a = 2x(15MHz), bandwidtns of 7 
MHz, interaction times of Ti = 117 /jsec., and chirp rates of 

/? = 2t( 3 x 10 10 ) rad/sec 2 . The spread spectrum code used for the experiment was a PN 
sequence of length equal to 63 chips. Furthermore, all 63 chips were contained within an 
information bit, thus providing a processing gain of 63. The chip rate of the PN sequence 
was 1.875 MHz. 

For a single tone interferer which was offset from the carrier by 180 kHz, and an 
interference—to-signal (I/S) ratio of 20 dB, the bit error rate (BER) for an energy per 
bit—to—noise power spectral density ratio (Eb/N 0 ) of 10 dB was 8 x lR 5 when the excision 
filter was present. Without the filter the BER was 0.5. The theoretical BER for BPSK 
systems without interference is 4 x 10- 6 . 


4.0 DIGITAL EXCISION TECHNIQUES 

In Section 2.0 a digital approach to interference removal w r as depicted in Fig. 

2(b). With respect to this approach, this section focuses on three areas: (a) the 
time-domain algorithms that can be used to calculate the filter parameters; (b) a review of 
the research which has been conducted in digital excision; and, (c) results obtained by the 
author on the problems of filtering stable and swept tone interferes. The signals in the 
sequel are all considered to be real. 

4.1 Preliminaries 

The time—domain algorithms referred to in Fig. 2(b) are based on time—modeling of 
stochastic processes. A well-known model [4],[27] is one in which the sampled signal x n is 
modelled as a p th order stationary autoregressive (AR) process defined as 

X n "h <UXn-l 4* — + apX ri -p = W n (6) 

where w n is w r hite Gaussian noise of variance aj. An AR process can be generated by 
passing samples of white Gaussian noise through a digital filter with a transfer function 
defined as 


H(z) =- l - - i— (7) 

1 + aiz* 1 4-.+ apZ-p A(z) 

where p is the order of the filter. 

For the moment, assume the objective is to whiten or, alternatively, to decorrelate a 
received set of samples obtained from a zero mean, stationary, p ttl order AR process. In 
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addition, assume that the second order statistics of the process are known. To whiten the 
signal, the set of parameters {ai, i=l,2, • • • ,p} must be determined, and a transversal filter 
with transfer function A(z) constructed with the a ; as tap weights. The data x n could then 
be passed through this filter whose output would be a white noise sequence. This is 
illustrated in Figs. 5(a) and (b). The parameters of the whitening filter A(z) can be 
determined in a straightforward manner by employing linear prediction concepts and the 
orthogonality principle [62]. 

The one step linear predictor is defined as 


Xn — "b &2^n-2 "h * ’ ’ " "4* &pXn-p- (8) 

The objective is to find those parameters a; which will minimize the mean squared error, 

E{f n 2 }, where f n is the forward prediction error x n — x n . The best estimate in a minimum 
mean squared error (mmse) sense is the one producing an error orthogonal to the data used 
in estimating x n ,i.e., 

EffnXn-i} = 0, i=l,2,- • -,p. (9) 


Substituting f n = x n — x n into Eq. (9) leads to the following p equations known as the 
Yule—Walker equations, i.e., 

Ti = airi-i + a2rj-2 4- • • • • + a,>Ti.p, i = 1,2, • • • ,p (10) 

where r, = E^nX^;} is the i th lag from the discrete auto correlation function of the AR 
process x n . This series of equations can be written in matrix form, i.e., 

Ra=r (11) 

where R is a pxp Toeplitz matrix whose elements are constant along the diagonals, i.e., 


Also 

and 



a = [ai a 2 • • • ■ apj* 
I = [rj r 3 • • • • r p ] 1 


( 12 ) 


(13) 

(14) 


where the superscript "t" refers to the transpose of a matrix. Equation (11) is knowm as 
the Wiener-Hopf equation. 


The solution to the Wiener-Hopf equation involves the inversion of the matrix R. 
It can be solved by several standard techniques, such as Gaussian reduction or Cholesky 
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decomposition. The number of mathematical operations required to obtain the inverse by 
these methods is p 3 /3 + ^(p 2 ) and p 3 /6 + ^(p 2 ), respectively. The term C{ ) refers to "in 
the order of ( ) operations, where an operation is defined as a multiplication or division 
plus an addition or subtraction. 

Alternatively, taking advantage of the Toeplitz structure of R, Eq. (11) can be 
solved most efficiently by means of the Levinson—Durbin (LD) algorithm [37] in ^(p 2 ) 
operations. This algorithm not only solves for the p th order equation, but also the lower 

order ones as well, i.e, for {i = 0,1,-,p—1). It is therefore an order recursive algorithm. 

It also solves for the minimum mean squared error, (defined also as <r 5 2 ) for the i th 
order. The orthogonality principle yields 

£min (i ) = E{fnX n } (15a) 



k=l 

= <Ti 2 (15c) 

The LD algorithm is listed in Table 1. It proceeds recursively to compute the 

parameter sets {ai ( %i 2 }, {ai( 2 )a 2 (% 2 2 },-,{ai(p)a 2 (p)-,a p (p)<r p 2 }. Note that a superscript 

has been added to the AR coefficients to denote the order. The final set 

{ai<p)a 2 (p)-,a p (pV p 2 } is the desired one and is used in the whitening filter in Fie. 5(b). 

The AR sequence x n can then be applied to the filter to yield a white sequence, f n , with 
variance <j 2 = <rj. 

P w 


The LD algorithm also implicitly calculates the parameters of a lattice filter 
[45],[46] as shown in Fig. 6. The coefficients of the lattice structure are the parameters 
{ai(»)= Ki , i=l,2,• • • ,pj. It can be shown that these parameters are less than one for a 
stationary process [61], and that this property immediately indicates that the roots of A(z) 
lie within the unit circle, which alternatively means H(z) in Eq. (7) is a stable filter. 

Table 1 

Levinson — Durbin Algorithm 


Initialization 


a ^ 1 ) = ti/tq 
(r 1 2 =(l-(a 1 d)) 2 )r 0 

Repeat for k=2.3. -,d 

{ k_1 

a k (k) = -^-! Ik “ ^ a i (k l)r k-l 

*k-l l i=l 

aiOO = ai(k-i) - a k (M ak.^-D, i=l,2, — ,k—1 

^k 2 = (1 - (a k < k )) 2 ) c k -i 2 . 
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FIGURE 6: WHITENING FILTER FOR A STATIONARY, pth ORDER AUTO 
REGRESSIVE PROCESS, IMPLEMENTED AS A LATTICE 
STRUCTURE. 
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The lattice filter is an interesting structure in that it contains the transfer functions 
of all lower order filters, i.e., Ai(z), A 2 (z),- • ■ *,A p (z), as well as the forward prediction 
error sequences {fi, n , i=l,2, • • • ,p} for each of these sections. For the p th order AR process 
under consideration, only the output sequence of the p th stage f^ n will correspond to the 
minimum mean squared error sequence. Another important pomt is that the lattice filter 
provides as a by-product a set of backward prediction errors denoted by {bi. n , i=l,2, - ■ - ,p} 
and defined as 


bi,n — X r 


-I 

k=l 


ak (i) Xn-i*k i=l,2,• * * ,p. 


(16) 


The backward prediction error is the difference between x..., and a weighted sum of i future 
samples {x n .i*i,x n .i, 2 ,- ■ • • ,x n }. For completeness, the forward prediction errors are written 
as 


fi,n = x n - Y ak^'Xn-k i=l,2,- - ,p. (17) 

k=l 

From the lattice structure, another definition for the forward and backward prediction 
errors is 

fin = f i- l.n Kibj-1 r>-1 (18a) 

bi.n = bj-i n-j — Kjfi-i.,). (18b) 

The backward prediction errors from each section of the lattice structure are orthogonal to 
one another. This is because the lattice can be viewed as carrying out a Gramm—Schmidt 
orthogonalization process [21] on the data x„. 

This latter property of the lattice can be appreciated by considering Eq. (18b). 
Assume bi-i.^i and fn. n , {i=1,2, ,p} are two vectors with an arbitrary angle 0 between 
them (in statistical terms, this implies they are correlated). The objective is to produce a 
new vector bi, n which is orthogonal to bj-i.^i by linearly combining bun-i and fn, n (in 
statistical terms, this implies the removal of the correlated portion of fi-i, n from bi-i.n-i). 

The parameter K, is the cosine of the angle between fi.],, and bj-i.^i and is equal to [29] 

E{ fi- l.nbj - 1 .n-l j 

Kj = - - (19) 

e{ f i- ln 2 j Ejbi-J.^,2 } 

Finally, the lattice parameters Kj are referred to as reflection coefficients because of 
the application of lattice filters to speech processing. The lattice is a model for the vocal 
tract, where each section of the lattice filter represents one tubular section of a certain 
cross sectional area of the tube [47]. As the acoustic wave propagates along the tract, 
portions of the wave are reflected at boundaries having different cross sectional areas. The 
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amount of reflection at these boundaries is measured by Ki. Thus, in analyzing speech 
waveforms, the Ki can be used to characterize a particular speaker. 

To this point, the solution to the Wiener—Hopf equation has been discussed for the 
case when the statistics of the process x n are known. In practice, these statistics are 
unknown. The next two sections discuss practical solutions. 

4.2 Block Least Squares Algorithms 

In the practical situation, only a sample realization of x n is available. The sequence 
of samples may have been generated by an AR process, or from any unknown process which 
had been assumed to be AR. Furthermore, only a finite set of samples are available. The 
problem of determining an optimum set of coefficients {a,, i = 1,2,- • • ,p} resides in the 
area of linear least squares theory [1],[34],[44|. These coefficients are obtained by solving a 
series of equations that arise from a mode which is similar to the AR process discussed 
earlier, i.e., each sample within a block of received samples {x n , n = 0,1,- - - *,T} is equal to 
a linear combination of the p most recent samples plus some error, i.e., 

x n = ajXn-i + a 2 x Ih2 4.+ apXn-p + f„. (20) 

The error f n is usually assumed to be white noise with Gaussian statistics. In Eq. (20), it 
has been assumed that the data has been fully windowed, i.e., x n = 0, n < 0 and 
n > T, as shown in Fig. 7. Pre—windowed and non—windowed cases are also shown. These 
latter two cases will be discussed later. 

Equation (20) is closely related to the AR process discussed earlier, except that the 
samples are deterministic values. Given the data, the objective is to determine the set 
{a i; i=l,2,* • • ‘,p} which will minimize the sum of squared errors between the estimate 
based on a weighted sum of past values, i.e., 

x r , = aiXn-i + ••■• + apx,,.,,, (21) 

and x n . For the fully windowed case, the following equations are obtained, written in 
matrix form: 



In matrix notation, this becomes 

e = x—Xa (23) 

where X has a Toeplitz structure. The optimum solution in a least sauares sense is the one 
producing a vector e orthogonal to the column space of X (see Ref. [if). In other words, if 
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FIGURE 7: THE EFFECTS OF THREE WINDOW TYPES WHICH CAN BE 
IMPOSED ON THE DATA. 
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the columns of X are considered as vectors in an N—dimensional Euclidean space (N = 
T+p-fl), and these p vectors form a sub—space within this N—dimensional space, then the 
optimum solution is that vector which minimizes the norm e'e. For this to be true, e must 
lie in the sub—space orthogonal to the sub-space created by the columns of X. Therefore, 
the optimum solution can be obtained from 

X*€=0. (24) 

Substituting for e in this equation, we obtain the following: 

X l [ Z - Xa ] = 0 (25a) 

PXa = X'x (25b) 

This is equ ; valent to Ra = r where R = X l X and r = X l x and are, respectively, estimates 

of R and r in Eqs. (It) and (14). For the fully windowed case, R is a pxp Toeplitz matrix, 
so that the LD algorithm can be used to solve for a. (Note, however that for an 

overdetermined set of N equations in p unknowns, forming R and r requires <?(Np 2 /2) 
operations, which can consume most of the computational effort for large N.) 

If the data are not windowed, so that no assumptions are made about the data 
outside the bloc! collected (i.e., the nonwindowed case, as shown in Fig. 7), there are 
several approaches to solving the problem, depending on thm assumptions that are made 
about the data (stationary or non-stationary), whether or not the roots of A(z) are 
required to lie within the unit circle, and general performance requirements. Three 
approaches will be discussed. 

A pproach 1 


This is the simplest of the three approaches. The problem involves minimizing the 
sum of squared forward prediction errors. The matrix equation of interest is 


r [p 1 


’ x p 


' x !>! •• 

X 

IV) 


[a,] 

* j>*i 

= 

X;,l 


Xp Xjr 1 •• 

■ X, 


a t 

. *t . 


. X T . 


Xji } 2 ’* 

X T-p . 


. a p. 


(26) 


Applying the same reasoning that was used in determining the optimum solution for Eq. 

(22), one obtains a matrix R which is no longer Toeplitz, so that the LD algorithm cannot 
be applied. Furthermore, the solution may not yield a filter with roots inside the unit 
circle. Although the efficient LD algorithm cannot be applied here, efficient algorithms 
nave been developed [28],[57] to provide a solution requiring t?(p 2 ) operations. 

One can derive a lattice structure (shown in Fig. 8) from the algorithm described in 
[57]. The reflection coeffic’ents in the forward and backward branches of each section (Kj f 
ana K ; b ) are not necessarilv equal, and their magnitudes are not necessarily less than 1. 
This latter property can be used to determine whether or not the roots of A(z) are within 
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FIGURE 8: LEAST SQUARES LATTICE STRUCTURE 
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the unit circle. As shown in [76], these parameters are bounded according to the relation 
0 < KifKjb < 1. Finally, the degree of disparity between Ki f and K, b is a measure of the 
amount of non—stationarity in the data [19]. 

Approach 2 (Burg’s Algorithm) [5].[60] 

The assumption is made that the process is stationary. Furthermore, it is required 
that the roots of A(z) be within the unit circle and that only the available data be used 
(i.e., the nonwindowed case). With these requirements, Burg developed an algorithm in 
which he suggested that the sum of squares of the forward and backward prediction errors 

be minimized with respect to the reflection coefficients {afb = Ki, i=1,2,-,p), i.e., 

minimize 


T 

^ ( fi.n 2 + bi. n 2 j, i = l,2,-,p. 

n=i 


(2T) 


The choice of Eq. (27) is based on the stationarity assumption for which the statistics of 
the forward and backward prediction errors are equal. Also, minimization of this equation 
with respect ot the i th stage reflection coefficients is equivalent to a local minimization of 
each stage of the lattice [45]. A local minimization of each section leads to a global 
minimization of the entire lattice. If one substitutes Eqs. (18a), (b) into Eq. (27), one 
obtains £, as a function of Taking its derivative with respect to and setting it equal 
to 0, yields the optimum solution for K,, i.e., 



i=l,2,- • -,p. 


(28) 


Equation (28) is the deterministic form of Eq. (19) for stationary processes. The Burg 
algorithm can be viewed as a solution to a constrained least squares problem, where the 
constraint is the imposition of the lattice structure in Fig. 6 on the data as opposed to the 
least squares structure in Fig. 8. This in turn forces a Toeplitz structure on the 
autocorrelation matrix. 


The initial conditions for the algorithm are fc. n = bc. n = x n . With this initial 
condition, Ki = af 1 ) can be calculated from Eq. (28). The block of data is passed through 
the first order lattice filter to produce the sequences f l k and bj.n, from which a^ 2 ) is 
determined, and so forth. After K p is calculated, one can determine the coefficients 
{afr), i=1,2, * •• - ,p} using the LD recursion if the power spectrum of the sequence x n is 
required. An analysis of the computational complexity of the Burg algorithm [29] indicates 
that 3Np — p 2 — 2N - p adds, 3Np - p 2 —N + 3p multiplications, and p divisions are 
required where N = T + 1. Storage of 3N + p + 2 values is also necessary. 
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Approach 3 


The constraint imposed on the data in the Burg algorithm has led to certain 
problems in its application to spectral estimation, such as spectral line splitting and biases 
in the frequency estimate [29],[48],[49]. These problems can be mitigated if £\ in Eq. (27) is 

minimized with respect to all coefficients {a k (p), k=1,2,-,p} [2],[48],[60],[82]. For the 

derivation of the solution, the forward and backward prediction errors are required: 


f f p 1 

fp.1 

_ 

HS • • 

- 

b 0 

bi 


- ^T-p 



A p 

x p *i 


Xo 

Xl 


Xf>l Xp-2 
Xp Xp-[ 


l T-p J 


X T-P*: X T-p-2 


XO 

Xl 


ai 

a 2 


X T-1 

X T-2 

X T-p J 1 

. a p J 

Xl 

X2 

4 . . V 

ai 

X 2 

X3 

• • • Xp.1 

a2 


l T J 


l P J 


(29a) 


(29b) 


These equations can be rewritten as 


§f = Xf — X f a 

(30a) 

§b = Xb 

(30b) 

The objective is to minimize 


£ p = + e h f e b 

(31) 


with respect to the set {a k 'p' 1 k= 1,2,* • • -,p}. This can be achieved by rewriting the 
prediction error equations as e = x — Xa. where 


§ = 

x = 



X = 



(32a) 

(32b) 

(32c) 


With these definitions, the solution is obtained from 

XtXa = X f x (33) 

An order recursive algorithm has been developed by Barrodale and Erickson [2], which 
utilizes the Cholesky decomposition method and certain shift invariant properties 
associated with the matrices X ( X and X f x The approach has good numerical qualities. 
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Marple [48] has also developed an efficient algorithm. The computational complexity of 
this latter algorithm for N = T + 1 is Np + 8p 2 + N + 7p - 8 additions, Np 4- 9p 2 + 2N + 
25p — 3 multiplications, and 5p +3 divisions. The algorithm needs N 4- 4p 4- 15 computer 
memory locations. 

Remark 


In approaches 1 and 3, it has been tacitly assumed that the inverse of X l X exists, 
i.e., it is of full rank. In many applications this will not be the case, particularly if there 
are fewer equations than unknowns (the under—determined case). For this case an infinite 
number of solutions for a exist which will minimize the mean squared error. Of this set of 
solutions, the one with minimum norm is typically selected. To obtain the minimum norm 
solution, the Moore—Penrose pseudo inverse [1] of XOC must be determined. 

4.3 Recursive Algorithms 

In the previous section a number of block structured algorithms 'were discussed 
regarding the solution to the Wiener-Hopf equation. These algorithms led to transversal 
or lattice filter structures. Furthermcie, with deterministic data, the problem was recast 
into a least squares one, in which the Wiener-Hopf equation was comprised of sample 
values of the process. 

In this section, recursive solutions to the Wiener-Hopf equation will be considered. 
In particular, stochastic gradient and recursive least squares algorithms for both the 
transversal and lattice structures will be discussed. 

4.3.1 Stochastic Gradient Algorithms 

The objective of the stochastic gradient approach is to make small adjustments to 
the filter coefficients in such a manner as to minimize the mean squared error between x n 

and x n . The adjustments to the filter coefficients are a function of the gradient of the mean 
squared error, which can be viewed as a multi—dimensional paraboloid with a unique 
minimum. This minimum point is sought by the algorithm by adjusting the coefficients in 
a direction opposite to the gradient. 

Gradient algorithms have been developed for transversal and lattice structures. 

Both structures will be discussed, with the former structure being used to present the basic 
concept behind the gradient technique. 

Transversal Structure 


If the statistics of the process are known, the recursion for updating the filter 
parameters is the following [84]: 

axu = ax - ^(VE{fpT 2 }) 


= a T + g(r - Ra T ) 
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(34a) 

(34b) 






where ax = [apt a 2 ,T ■ • • a^x] 1 , T is the time index, ^ is a small positive constant which 
determines the rate of convergence of ax to its optimum value a op t, and V is the gradient 
operator, i.e., 


d/da, 

d/da 2 


L d ! d \ 


(35) 


The constant /x must be selected within certain bounds to ensure convergence of the 
algorithm to the optimum solution. The bound is defined by 0 < /x < 2/^ max where A max is 
the maximum eigenvalue of the auto correlation matrix of the process. 

Practical implementation of the algorithm involves approximating r and R in Eq. 
(34) by lx and Rx-i at time T, where 


XT-1 
XT-2 


x tXt-1> 


(36) 


XT-P 


and Rx-i = xt-iXt-i 1 - After substituting r T and Rx-i into Eq. (34b), the update equation for 
the taps of the transversal filter becomes 


where 


§t.i = in + /xfp tXt-i 


(37) 


fp,T = x T - ax^x-i. (38) 

It should be noted that the term f p .xXT-i is a noisy, unbiased estimate of the gradient 
E{fp iT xx-i}. Eq. (37) is often referred to as the Widrow-Hoff LMS algorithm. 


Remarks 


(1) With the above approximation, the gradient is noisy; however, E{ax) —♦ a op t if the 
condition 0 < fi < 2/\ wax is met. 

(2) The tap weights fluctuate about a op t- The degree of fluctuation is a function of the step 
size /x. Small /x implies averaging over more past data and therefore the variance of the 
fluctuations is smaller. The cost is a slower convergence rate. 

(3) Tap weight fluctuations lead to an excess error in steady state. Thus the mean squared 
level of the output error sequence f p x is larger than the minimum achievable. This leads to 
an excess mean squared error. 
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(4) The average mean squared error converges to its steady state value if and only if the 
step size parameter /i satisfies the condition 

0 < < - -= —(39) 

^ p pr 0 ' ' 

X* 

i=l 

for a stationary process where ro is the power in the signal. If this condition is met, the 
LMS algorithm is convergent in the mean square. This imposes a tighter bound on /*. 

(5) In certain applications such as speech processing, where power levels can fluctuate 
significantly, the LMS algorithm can become unstable. In these cases the step size is 
usually made a function of time [23]. One possibility is [22] 

H = —r—• (40) 

pf 0.T 

In Eq. (40) 0 < a < 2 and fo.T is an unbiased estimate of the power in the sequence xt, i.e., 
? o.t = (!—/?)r 0T .i + /?x t 2 (41) 

where 0 < 0 < 1 and fo.-i is an estimate of the initial condition for ro at T = -1. 

(6) The rate at which the LMS algorithm converges to the steady state value is dependent 
on the eigenvalue spread. It has been determined [16],[20] that the overall time constant r a 
for any tap weight a, is bounded by the relation 



(7) The LMS algorithm can track processes which are non—stationary. The degree of 
success in its tracking capability is a function of the amount of non—stationarity in the 
process, the eigenvalue disparity and the step size fi. 

Lattice Structure [17],[16],[23] 

The derivation of the gradient adaptive lattice algorithms is similar to the method 
used to derive the LMS algorithm, except that the gradient is applied to each reflection 
coefficient. There are two approaches which can be taken. If the data are assumed to be 
nonstationary, then the forward and backward reflection coefficients K, f and Kj b are chosen 
to minimize, respectively, the mean squared errors E{f s t 2 } and E{bi,t 2 } of each lattice 
section. The statistical update equations are 

K i T *i f = K iT f - 2 /Zj f V f E{b i t 2 } 

K,T.i b = K,t*> - j gi b V'>E{fi x 2 } 
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(43a) 

(43b) 






where gi f and n i b are the step sizes for the i th section, and V f and V b are the 
one-dimensional gradient operators, d/dK J and d/dK\ h , respectively. After substituting 
for fij and bi.Tfrom Eqs. (18a),(b) into Eqs. (43a),(b), taking the derivative and 
approximating the expectation operator by the instantaneous values, the following 
approximate gradient equations are obtained: 

Ki, T *i f = K,.t* + /*i f fi-i,Tbi.T (44a) 

Ki,T*i b = Kj,T b 4- fi\ b fi.Tbi-i.T-i- (44b) 

The second approach assumes the data are stationary, implying that Ki f = Kj b = 
Ki. The objective is to find that K; which minimizes the sum Ejfi.T 2 } + E{bi,r 2 }. The 
approximate gradient equation is 

K lT u = Ki.x + Hi[ fubi-i.T-i + fi-abi,T ]• (45) 

Remarks 


(1) The gradient lattice algorithms converge to the steady state faster than the LMS 
transversal algorithm [12],[17],[23] for the case in which there is a large spread in the 
eigenvalues. This is because the lattice algorithm has more degrees of freedom in adjusting 
its parameters; it has at least one convergence parameter g for each section. This should 
be compared to the LMS algorithm, in which there is only one parameter for the entire 
filter. 

(2) The step size can be power normalized, in a manner similar to that for the LMS 
algorithm [17],[23]. For example, in Eq. (45) normalization is of the form 


RT = l/°'i.T 2 


(46) 


where 


^ i.T 2 = (1— f3)h,T-l 2 + /?[ fi. T 2 + bj.T-! 2 ] (47) 

and 0 < 0 < 1. 

(3) The gradient lattice algorithm is considerably more difficult to analyze than the 
gradient transversal algorithm because of the highly nonlinear nature of the adaptation 
(i.e., each section converges at a different rate). It has been found empirically [23] that for 
the power—normalized case and the stationary form of the lattice structure, stabilization 

will be maintained if 0 < l/c’; 2 < 1. 

4.3.2 Recursive Least Squares Algorithms 

In this section the pre—windowed (Fig. 7) recursive least squares algorithm [7] will 
be discussed. The non—windowed case has a recursive solution as well, albeit the algorithm 
is more complex and requires a more sophisticated start—up procedure 
[12],[13],[35],[41],[64], 
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The matrix equation for the pre—windowed case at time T is 


o —* 

> -— 

= 

'V 

X 1 

— 

' 0 
x o 

0 

0 

• • • 0 ' 
... o 


Fa i ' 

a 2 

(48) 



. X T . 


. X T-1 

X T-2 

• • • Y 

x T-p . 


a 

L Q P 



The optimum solution for a p j for a block of data where T > p can be solved using, for 
example, the generalized Levinson recursions in [57], modified for pre—windowed signals. 

The major computational burden associated with block methods of computing the 

optimum coefficients is the matrix R = XOC, as noted earlier, when the number of samples 

is large. Furthermore, R must be recomputed when the least squares estimate is to be 
updated when a new data sample arrives, i.e., as T increases. A number of algorithms over 
the years have evolved which address these problems. 

The first algorithms were recursive in time for a specific filter order. The original 
one [15], based on a Kalman filter model of the least squares problem, required ^(p 2 ) 
operations per time update. A more efficient method [11],[36] known as the fast Kalman 
algorithm was developed requiring t?(8p) operations per update for AR modeling. The 
number of operations was reduced to <?(5p) by the FAEST algorithm [6],[7] i.e., fast a 
posteriori error sequential technique, and FTF algorithm [8], i.e., Fast Transversal Filters 
algorithm. The fast Kalman and FAEST algorithms will now be discussed. 

The optimum solution at time T for a p ,h order system, derived from Eq. (48), is 
obtained from R t1 It = It where 


T-l 

Et-i = ^ XiXi 1 (49a) 

i=0 

T 

lT=^x i Xi. 1 t (49b) 

i=0 

Si = [xj Xu-Xi-p.jp. (49c) 

The problem is to determine a™ when the sample xj.j arrives. The key to the development 

of a sequential algorithm is the existence of simple time updating formulas for Rj and r t*i 

in terms of Rt-i and r T , respectively, and the most recent data {xt.i,*t}- These update 
formulas arise from the Toeplitz structure of the data matrix X in Eq. (48). Using these 
updating formulas [57], it can be shown that the linear prediction update equation is of the 
form 


It.] = It + Rt '’xTfp.T.j 


(50a) 
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= ax + Kg.t^.t*!- 


(50b) 


where Kq.x is referred to as the Kalman gain vector and f£,x,i is the a priori error defined as 

fgx i = x T »! - a^xx- (51) 

It is sometimes referred to as the predicted error since the prediction of xx*i, namely ax^x, 
is based on the tap weight vector at time T, not T+l. The fast Kalman algorithm [11], 
[36] is a series of operations per time step concerned with updating the Kalman gain vector 
Kg,t- 

With some modifications to the fast Kalman algorithm, the FAEST algorithm [6], 
[7] was developed using an alternative formulation, which had the effect of reducing the 
number of operations required to update the Kalman gain vector. The update equation in 
this case is 


ax*i = ax + Rx-f I XTfp,T*i 


(52) 


= ax + K(xxfp.x*i 

/ 

where Kg.x is the alternative Kalman gain and 

fpT*i = x T *i - ax.^xx (53) 

is the a posteriori error. 

When either of these algorithms is started at T = 0, the rank of the autocorrelation 

matrix Rx is less than the filter order p. To prevent a singularity from occurring in Ry for 
T < p, the initial minimum mean squared error is set equal to a small constant 6. This is 
equivalent to initializing the algorithm’s autocorrelation matrix to <51, where I is the 
identity matrix. This initialization results in a slightly biased least squares solution. 

Remarks 


(1) The inverse matrices in Eqs. (50) and (52) have the effect of equalizing the eigenvalues. 
Thus, these fast least squares algorithms are independent of eigenvalue spread in contrast 
to the LMS algorithm. The smaller the eigenvalue spread, the smaller the difference 
between the LMS and fast algorithms. 


(2) If the right hand side of Eq. (48) is premultiplied by A 1 / 2 where 


' A T 0 0 

0 A 1 ' 1 - • • 0 


0 0 • • • 1 


(54) 
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and A < 1, then one is concerned with a weighted least squares problem in which the old 
data is weighted out. This provides the capability of tracking time-varying processes. 

The memory of the algorithm in this case is approximately 1/(1—A) samples. Furthermore, 
the initial starting condition of £1 becomes insignificant for large T. 

(3) The fast algorithms have a problem with numerical stability due to roundoff error 
accumulation [8],[40]—[41]. Rescue techniques have been developed [8],[41] to deal with this 
problem. Cioffi [8] also developed normalized algorithms which increase the dynamic range 
performance of finite precision algorithms, leading to a reduction in the number of rescues 
that are required. The choice of the initial value 6 can also lead to instability [8]. To 
remove the requirement for this initial condition, Cioffi developed a fast start-up 
procedure for the FTF algorithm based on the minimum—norm least squares solution [1], 
This procedure produces the exact least squares solution mom T = 0. 

(4) Both unnc r malized and normalized recursive least squares algorithms have been 
developed for the lattice structures [12],[13],[35|,[43],[56j,[58],[59J. Though more stable 
than their transversal counterpart, they can still develop numerical problems [8] and 
eventually diverge after continuous operation. One advantage of the lattice structure is the 
availability of all filter orders up to some upper limit. This may be useful in some 
applications in which the correct order to use is unknown. The disadvantage is the increase 
in the number of computations/iteration. Table 2 summarizes the computational effort per 
iteration required by several of the recursive algorithms [8]. 

Table 2 

Computational Effort for Linear Prediction 



Algorithm 

Operations 

LMS 

2p 

Fast Kalman 

8p, 2 div. 

FAEST 

5p, 2 div. 

FTF 

5p, 2 div. 

Normalized FTF 

9p, 2 sq. roots, 2 div. 

Unnormalized 


Lattice 

lip, 3p div. 

Normalized Lattice 

18p, 3p sq. roots, 


4p div. 


4.4 A Review of Research Work in Digital Excision 

In the previous section, several adaptive algorithms and their respective 
characteristics were discussed. With the advent of high-speed digital technology, these 
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algorithms are becoming feasible in interference excisors. Many of these algorithms have 
been evaluated by researchers either analytically, under computer simulation, or through 
implementation. 

The earliest work revealed that Wiener filter theory could indeed be applied to the 
suppression of interference in DS/SS systems, thus providing significant improvement in 
SNR. This revelation spawned a great deal of research in three areas: (a) algorithm 
performance under various interference conditions; (b) the development of theoretical 
performance models; and (c) filter structures (linear interpolation error filters and filters 
employing decision-feedback). This section reviews a large portion of this work in 
somewhat chronological order starting with the basic principle. 

The filtering of narrowband interference from DS/SS systems is based on the 
property that the signal and receiver noise are noncoherent and therefore unpredictable, 
whereas the interference is coherent over the short term and hence predictable. With the 
linear prediction algorithms of Sections 4.1 to 4.3, this coherent component can be 
estimated and removed from x n in Fig. 2(b) using transversal or lattice filters. The output 
of the filter is an error signal f n consisting of the sum of the filtered versions of the desired 
spread spectrum signal, thermal noise and residual interference. 

Hsu and Giordano [24] appear to have been the first to suggest this method of 
interference removal. Their input signal model was a sequence x n = s n 4- n n + i n where s n 
was a chip sample of amplitude ±1, n n was a sample of white Gaussian noise of variance 
o nse 2 and i n was a sample of narrowband interference with autocorrelation function p(i), i = 
0, ±1, ±2, • • •. The assumption was made that the signal, noise and interference were 
independent from one another and the signal sequence was white with autocorrelation 
function 


R ss = <5(n) = 1, for n = 0, (55) 

= 0, for n i 0. 

With these assumptions, the interference becomes the predictable component as now 
shown. 


Consider a p th order linear predictor as in Eq. (8). Its coefficients are determined 
from Eq. (11). The elements of the autocorrelation matrix, R, are given by 


r 0 = E{x n 2 ) (56a) 

= 1 + <T nse 2 + p{ 0) (56b) 

where p( 0) is the interference power, and 

r j = E{xnx rhi } (57a) 

= Ejinin-i} (57b) 

= P(i), i = ± 1, ± 2, • • •, ±p. (57c) 


Equation (57b) reveals that the interference is the predictable component, since the 
autocorrelation function r* depends only on /?(i). 

The output of the excision filter is the error sequence f n defined as 


- 27 - 




(58a) 
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(58b) 


The first two terms are non-white signal and noise sequences, and the third ’S the residual 
interference. The direct sequence despreader correlates the error signal with the PN 
sequence, as shown in Fig. 2(b). 

Hsu and Giordano studied the performance of prediction (i.e.,single—sided) 
transversal filters with respect to the suppression of interference modelled as the sum of 
sinusoids of equal amplitude and random phases. Furthermore, the interference spanned 
20% of the chip rate. In their analysis, the filter coefficients were calculated by means of 
the LD and Burg algorithms. 

The receiver’s performance was measured in terms of a signal—to—noise (SNR) ratio 
improvement parameter at the output of the PN correlator, i.e., the ratio of the SNR with 
whitening, to the SNR without whitening. If the output of the bit detector is the test 
statistic U, then the SNR ratio when the excision filter is present can be shown to be equal 
to 


SNR 0 = 


E 2 {U} 
var { U} 


(59a) 


L 


L 2 


P 

p p p 1 


X aj2 + 

/?( o ) - 2 y aj ,(j) + y y aja ^(j-k) 

mmd tom mb 

+ 

j=l 

L J=1 j=l k=l J 




P 'i 

iTn.se 2 aj 2 

J=! ' 
(59b) 

(59c) 


where L is the number of chips/bit, <r 2 is the overall noise power and the chip rate is 1 Hz. 
(It should be noted that for chips of amplitude v c and duration T c , the energy per chip E c 
is v^c. The energy per bit E b is therefore LE C . For chips of unit energy, E{U} = L.) 
This noise power consists of three parts, as noted, in the denominator of Eq. 59 (b): a 
signal distortion term due to the excision filter, residual interference, and Gaussian noise. 


Without the filter, the SNR is simply 


SNR no = 



(60) 
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Thus, the figure of merit for the excisor is the ratio 
V = SNR 0 /SNR no , i.e., 


V = 


L(/?(0) + Cnse 2 ) 
a 2 


(61a) 


Mo) 

~-, for p(0) >> <Tn.se 2 - (61b) 

<T 2 

Both the LD and Burg algorithms yielded substantial SNR improvement for large 
interference-to—signal (I/S) ratios. For example, for L=10, <r.^ e 2 — 01, 10 tones and SNR no 
= —16 dB (large p{ 0)), T) was approximately 20 dB for the Burg algorithm and 12 dB for 
the LD algorithm. The LD algorithm demonstrated a leveling off in 7j for large I/S, 
whereas the Burg algorithm did not. The inferior performance exhibited by the LD 
algorithm can be attributed to its characteristic tendency of providing a larger bandwidth 
notch compared to the non—windowed least squares algorithms, resulting in more residual 
interference. Finally, performance did not vary much with filter order for nigh I/S 
conditions, i.e., low ( p = 4) and high order (p = 29) filters produced very similar results 
for the case of 10 tones. For 100 tones the difference started to manifest itself, with the 
high order filter (p = 29) providing better performance, suggesting the removal of more 
interference. In contrast, for low I/S conditions, the performance of the high order filter 
was worse than the low order one. This difference in performance was due to the 
additional signal distortion intro^ ’ced by tne extra zeroes in the larger order whitening 
filter. 


Li and Milstein [38] compared the performance of the single-sided (prediction error) 
filter to the two-sided (interpolation error) linear phase transversal filter with the same 
number of taps as the prediction filter. The two-sided filter configuration is shown in Fig. 
9. Single and multiple tone interference were considered. 

For the case of a single tone interferer, an expression for SNR improvement was 
derived using the filter coefficients obtained from the exact solution to the Wiener—Hopf 
equation. The results showed that when the tone was in close proximity to the carrier 
(approximately within .05 Hz for a normalized sampling rate of 1 Hz) and the noise 
variance was zero (the ideal case), the two-sided filter provided up to 6 dB additional SNR 
improvement over the prediction filter. This was due to the more accura-e estimate of the 
interference produced by the interpolator. For frequencies further away from the carrier, 
the performance difference was not significant, suggesting the interference estimates 
produced by the interpolator filter were less accurate, likely due to the decrease in 
correlation between interference samples for increasing frequency offsets. 

Analytical expressions were not developed for the multi- tone case because of the 
difficult nature of the problem. However, simulations were conducted with the 
Widrow—Hoff LMS algorithm. In general, the two-sided filter performed better than the 
prediction filter by 3 to 10 dB for I/S’s ranging from 2 to 30 dB. 
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Li and Milstein [39] also analyzed the performance of a direct sequence QPSK 
spread—spectrum system using complex prediction and interpolation error filters in the 
presence of pulsed CW interference. They examined the convergence behaviour of the LMS 
algorithm during the "on" and "off" portions of the interference, and found that for the 
case when the product of the number of taps times the jammer power was much greater 
than the sum of the signal and noise powers, i.e., p/?(0) >> S + 0 - nse 2 , the taps changed 
much slower during the "off" portion of the pulse, thereby maintaining good rejection 
capability when the "on" portion of the interference returned. Furthermore, by deriving an 
autocorrelation function for the pulsed signal, they were able to calculate the optimum 
weights and, therefore, obtain an expression for SNR improvement for the worst case 
situation, i.e, when the start of the pulsed interference coincides with the leading edge of a 
symbol. This situation generates large interference spikes at the output of the filter, which 
gradually diminish as the interference cancellation begins to take effect. For both filter 
types, 7] demonstrated a large degree of sensitivity to filter order for a small duty factor 
(1%), and almost no sensitivity for larger duty factors (10%). This behaviour suggests a 
requirement for higher filter orders to attenuate interferes of increasing bandwidths. 

Ketchum and Proakis [32] based their work on that of Hsu and Giordano. They 
considered: (a) single—band and multiple—band interference; and, (b) two additional 
algorithms for estimating the filter coefficients. The interference consisted of either 
narrowband Gaussian noise or a sum of closely spaced sinusoids of equal amplitude and 
random phase, spanning 20% of the chip rate in a single-band segment or in multiple-band 
segments. Besides the LD and Burg algorithms, they also assessed the performance of the 
non—windowed unconstrained least squares algorithm [Eq. (33)] and a frequency domain 
algorithm based on the FFT. The next two paragraphs provide a brief description of this 
latter algorithm. 

The frequency domain method relies on the characteristic that the power spectral 
density of the PN sequence is relatively flat, while the spectrum of narrowband interference 
is highly peaked. The first step in this method is to estimate the power spectral density of 
the received signal. The spectral estimate can be obtained, for example, by the Welch 
algorithm [65]. Once the power spectral density of the received signal is estimated, the 
interference suppression filter can be designed. 

The transfer function of the suppression filter for a p-tap filter is defined as the 
reciprocal of the square root of the power spectral density at equally spaced frequencies 
(0, R. c /p, 2R c /p,- • • •, (p—l)R c /p), where R n is the chip rate. This transfer function is 
defined as 


A(n) =- - - ■ —. e J r [ 2 , n = 0,1,-• • ,p-l (62) 



— fPlhn 

where P{ } is the power spectral density function. The term e " J p ^ 2 / n produces a 
realizable filter with delay (p - l)/2 samples. Using the inverse DFT on A(n), the impulse 
response (tap weights) of the filter can be calculated. The result is a linear phase finite 
impulse response (FIR) filter, which yields an interpolation error filter as discussed earlier. 
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Theoretical and simulated performance curves were presented for both the time 
domain and frequency domain algorithms. 

For the time domain algorithms, SNR improvement for single—band interference of 
a specified bandwidth was very similar for 4 th order and 15 th order filters, suggesting that 
accurate modeling of the interference was not critical. In contrast, for multiple band 
interference, the filter order that was required depended on the number of interfering 
bands; at a minimum, one pair of zeroes was recommended per band for BPSK signaling 
(or one complex zero for QPSK signaling). Finally, the suppression filters exhibited 
little difference in performance between the two models representing the interference, i.e., 
whether it was modelled as narrowband Gaussian noise or as a sum of sinusoids. 

They also considered the prediction error filter A(z) in cascade with its matched 
filter A*(l/z’), as shown in Fig. 10. This arrangement produces a linear phase filter like the 
interpolation error filter, but whose center tap is not equal to unity. This structure also 
resembles a maximum likelihood (ML) receiver (to be discussed later), in which the 
received signal is whitened by the error filter A(z) and the result is applied to its matched 
filter. As an example of its performance, this structure exhibited approximately 5 dB of 
additional improvement in performance over its prediction filter counterpart for the case of 
I/S = 20 dB. 

The frequency domain algorithm also proved to be a viable means of interference 
suppression. However, there were two main differences between it and the time—modeling 
approaches. First, the notch of the frequency domain filter was not as deep, leading to a 
degradation in tj of a few dB for the example considered. Second, the time domain 
modeling techniques required much fewer data samples than the frequency domain 
approach to provide good interference suppression. For example, the former provided 
suppression approaching that predicted by theory with as few as 50 data samples, whereas 
the Welch algorithm required 992 points to obtain similar performance. This latter result 
may prove to be unacceptable in an interference environment which changes quickly. 

As a final point, Ketchum et al. derived BER expressions based on the assumption 
that the noise at the output of the adaptive filter was Gaussian. (Recall from Eq. (59b) 
that the noise consisted of three terms, two of which were non—Gaussian.) The resulting 
BER curves compared favourably to those obtained from a Monte Carlo simulation, 
confirming the validity of the assumption. Furthermore, considerable improvement in 
BER performance could be achieved, as noted earlier, by passing the output of the 
prediction error filter through its matched filter. For a processing gain of 60, 

I/S = 20 dB, 10 sinusoids spanning 209c of the chip rate, and E n /N 0 = 10 dB, the BER 
improved from 2x10-" (unmatched) to IxlCH (matched). 

Saulnier et al. [70]—[74] emphasized the implementation of the LMS algorithm as 
applied to interference suppression. They considered both charge transfer devices (CTD’s) 
and digital hardware. Since direct implementation of the LMS algorithm requires two 
multipliers per filter tap (i.e., one for the tap weight update increment, and a second to 
calculate the signal sample), a burst processing structure was used in order to reduce the 
hardware multiplier count. This approach allows implementation of the adaptive 
algorithm with only two multipliers regardless of filter order. The cost, however, is a 
reduction in bandwidth processing capability by a factor equal to the number of taps. 
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FIGURE 10: WHITENING FILTER FOLLOWED BY ITS MATCHED FILTER. 
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In the CCD implementation [70], a stable tone jammer was considered for cases 
where I/S = 10 and 15 dB. A 7 chip PN sequence with a 26.3 kHz chip rate (sample rate) 
was used along with a 16 tap interpolation filter (8 taps on each side of the center tap). 
Empirical BER curves were compared with a theoretical curve based on SNR 0 at the 
output of the PN correlator. The results wpre interesting, in that for Eb/N 0 = 10 dB, the 
theoretical result yielded a BER of 2xl0- 5 for both interference conditions, whereas the 
experimental results were 2xl0- 3 (I/S=10 dB) and lxlO- 2 (I/S = 15 dB). The main reason 
for the discrepancy between the theoretical and experimental results appears to be the use 
of a 7 chip PN sequence relative to the number of taps in the filter. This would lead to a 
certain degree of code predictability for cases in which adjacent bits were equal, resulting in 
some code cancellation. Other reasons for the discrepancy would be due to (a) the use of 
the LMS algorithm, which only approximates the Wiener—Hopf equation, (b) inaccuracies 
due to implementation, and (c) charge transfer inefficiencies in the device. 

A digital implementation of the LMS algorithm using the same burst configuration 
was assessed [71],[73]. The sampling rate in this case was increased to 177 kHz. In 
addition, the input was quantized to 8 bits and both 7 chip and 31 chip PN sequences were 
used; the full PN code was contained in one symbol. Three interferer types were examined, 
i.e., a 30 kHz stable tone, narrowband Gaussian noise of bandwidths of 1 and 2 kHz 
centered at 30 kHz, and a swept tone [73]. 

With respect to the stable interferers, the effect of interference bandwidth was quite 
evident. For the case of I/S=15 dB and E b /N o =10 dB, the tone jammer yielded a BE?, of 
2xl0 4 ; the 1 and 2 kHz noise jammers yielded BER’s of lxlO- 3 and 5xl0‘ 3 , respectively. 

This performance degradation stems from two phenomena: (a) the wider notches created by 
the excision filter, which has the effect of introducing more signal distortion, and; (b) a 
reduction in the predictability of the interference as the interference bandwidth increases. 

In addition, for the single tone case, the effect of the level of interference power for fixed 
signal power with and without the excision filter was examined. Without the filter, and 
i/S ranging from 0 dB to 20 dB, the BER ranged from 2xl0 4 to 0.5. With the filter 
present, BER ranged from lxlO 4 to 5xl0 4 , indicating the sensitivity of the system without 
the filter, to different power levels. 

With respect to the swept tone interferer, the objective was to assess the effect of p 
in Eq. (37) for several sweep rates. The sweep rates ranged from 0 to approximately .001 
Hz/sec. (normalized with respect to the sampling rate of 177 kHz), and fi ranged from 0.002 
to 0.029. For small fi (0.002) and a sweep rate of .001 Hz/sec., the BER was 0.15; 
decreasing the sweep rate to 0.0001 Hz/sec., resulted in a slightly less BER, even though 
the sweep rate was reduced by a factor of 10. On the other hand, for larger g (0.029) and 
the same sweep rates, the BER’s were 5xl0 2 and 8xl0' 3 , respectively, confirming the 
requirement for a larger convergence parameter in the LMS algorithm in order to better 
track the tone. 

litis and Milstein [25] developed BER expressions in terms of the excision filter 
coefficients for two receiver types. The first receiver consisted of a bandpass RF section 
followed by an analogue tapped delay line; the second has already been discussed (Fig. 

2(b)). Each receiver was subjected to either a tone or narrowband Gaussian jammer and 
evaluated in terms of its BER performance. Three design criteria were used to calculate 
the filter coefficients. 
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The first criterion used an autocorrelation function derived from the combination of 
signal, noise and jammer. This refers to the familiar case already discussed in which the 
least squares filter whitens the spread spectrum signal, jammer and noise in a proportional 
sense, with most of the whitening being applied to the jammer. 

The second criterion used only the nc ; ce and jammer in the autocorrelation function 
(it was assumed that somehow the signal could be removed). In fact, the only difference 
between the autocorrelation matrices of criteria 1 and 2 was the disappearance of the signal 
power from the diagonal elements of R in Eq. (56b) when criterion 2 was applied. The 
second criterion was motivated by the structure of the Maximum Likelihood (ML) receiver 
[83] for the problem of the optimum detection of a signal in coloured noise n n c as shown in 
Fig. 11. The ML receiver consists of a whitening filter for the coloured noise (in this case 
noise plus jammer) followed by a filter matched to the whitening filter and the signal. 

The third criterion used tap weights which would produce zeroes of A(z) located on 
the unit circle. This is a limiting condition of the analytical equation developed by Li and 
Milstein [38] for a stable tone jammer. As the jammer power approaches infinity, the 
zeroes of the transversal filter tend towards the unit circle. 

The BER results produced by these three criteria were compared to each other for 
both prediction and interpolation filters. They were also compared to the performance of 
ML receivers derived for tone and Gaussian noise jammers. The ML receivers provided a 
benchmark for the other receivers. Of the three criteria, the second criterion, in which the 
signal power did not form part of the autocorrelation matrix, yielded the best results, 
approaching the performance of the ML receiver. The filter consisted of the prediction 
filter followed by its matched filter. This structure is similar to the one proposed by 
Ketchum in Fig. 10, in which the filter taps were calculated using the first criterion. 
Achieving the performance level of the ML receiver would be difficult in practice unless the 
signal could be removed from the input x„. 

litis suggested one such approach [25], employing decision—feedback, as shown in 
Fig. 12. The spread signal contained in x,„ delayed by one bit (L chips), is removed most 
of the time for low BER’s. The filter coefficients are thus calculated on the basis of the 
jammer and noise. The basic concept will now be described. 

Consider a linear interpolator and stationary processes. The interference i^L is 
estimated from "p" future and "p" past samples from the sequence yn-L- The mean squared 

value of the error (call this error f^i/) between y, h L and its estimate i„-L is minimized with 

respect to the set of 2p coefficients {a i; i = -p,-1,1,- • • ,p}. The optimum set of 

coefficients is then used in the excisor to which x n is applied. Since it was the mean 
squared value of the error sequence fnV which was minimized with respect to the 
coefficients a, and not f n , the signal s n in the sequence f n will not be as distorted by the 
filtering operation. 

In another paper by Ketchum [31], three decision feedback schemes were suggested. 
The basic principle is illustrated in Fig. 13(a). The concept is somewhat different from 

litis’. A delayed chip decision s„-i is fed back and subtracted from a delayed version of the 
input x n , producing the sequence yn-i = ip-i + n th |. This sequence is applied to a linear 
prediction filter G(z) = 1 - A(z) (A(z) is the error filter as defined in Eq. (7)) which forms 
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FIGURE 11: MAXIMUM LIKELIHOOD RECEIVER FOR THE CASE OF SIGNAL 
IN COLOURED NOISE. 
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FIGURE 12: 


r n - L 


EXCISION FILTER EMPLOYING DECISION FEEDBACK, WHICH 
APPROXIMATES THE MAXIMUM LIKELIHOOD RECEIVER. 


- 37 - 







an estimate of the interference i n . The optimum estimate is obtained from the 

orthogonality of the error f n (which is the difference between x n and i n ) to the data set 
{yn-i, i = 1.2,* • • ,p}. The performance of this method is sensitive to the chip error rate. 
Under the condition of no chip errors, the signal distortion term in Eq. (59b) is zero. For 
low signal—to—noise—plus—interference conditions, the performance tends to drop off 
quickly because of the increase in chip errors. Two other techniques, which provide better 
performance at the expense of system complexity, are described next. 

An improvement to the scheme in Fig. 13(a) is shown in Fig. 13(b), in which an 

evolving estimate for the j th bit djk is quantized at each chip interval kT c and multiplied 
by the local PN correlator, yielding improved estimates of the signal. This technique 
produces more chip errors near the beginning of the bit interval, but as the bit evolves at 
the output of the summer, the chip error rate decreases. The improvement in SNR that 
one achieves for this case thus becomes a function of the processing gain. The performance 
of this technique rivals that produced by the whitening filter/matched filter combination 
discussed earlier (Fig. 10). 

Better performance is achieved with the more complex configuration shown in Fig. 
13(c). The sequence y^i, with the spread spectrum signal having been "removed" by the 
technique in Fig. 13(b), is applied to an interpolation filter R(z). This filter is formed by 
first combining the error filter A(z) = 1 + G(z) with its matched filter A’(l/z‘) and then 
removing the center tap of value g'(0). The output of R(z) is the interpolated estimate of 

the interference in-j>i' which is removed from a weighted and delayed version of the input 
g / ( 0 )x n . p .i. As an example of the performance of this approach, for I/S = 20 dB, an 
interferer with bandwidth equal to 20 % of the chip rate, and a processing gain of 60, this 
system experienced increases in 77 of 14 and 7 dB over the systems in Figs. 10 and 13(b), 
respectively. 

Other research has been conducted on decision feedback approaches. For more on 
this topic, one can refer to the work of Takawira and Milstein [80] and Shah and Saulnier 
[75], 

To this stage the emphasis in adaptive digital excision techniques has been on 
transversal filter structures. Gradient lattice [19],[74] and least squares lattice [ 66 ],[33] 
structures have also been examined. In [19], the gradient lattice was used in filtering a 
stable tone; in [74] it was implemented on a TMS32010 processor. 

In [ 66 ] several receiver configurations employing least-squares (LS) lattice filters 
were assessed with respect to excision and channel equalization. The interference was 
coloured noise. Two important issues were addressed: (a) utilizing the LS lattice forward 
predictor in series with the backward predictor (prediction filter/matched filter 
combination) to remove the interference; and, (b) demonstrating that the LS lattice could 
be implemented as a joint process estimator to carry out decision—directed equalization and 
interference suppression simultaneously, provided the SNR was low (which can be the case 
in DS/SS systems) and the signal spectrum was reasonably flat. It was shown that the 
joint process estimator approximated a transfer function consisting of two filters in series, 
i.e., a filter matched to the channel’s transfer function and a filter for whitening the 
coloured noise. When the joint processor estimator was used, simulation results for \ = i.O 
in Eq. (54) produced results which were marginally better than the LS lattice whitening 
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filter, even when the LS whitening filter was preceded by the channel’s exact matched 
filter. However, for A < 1.0, the use of simultaneous equalization/whitening provided 
inferior results compared to the LS whitening filter. This result suggests that the joint 
process estimator is more sensitive to the misadjustment error for cases when A < 1.0 (see 
Section 4.5.2 for the definition of misadjustment error as applied to the recursive least 
squares algorithm). 

4.5 Adaptive Filtering of Stable and Swept Tones 

This section presents theoretical and simulation results of the application of least 
squares filters to the removal of interference from a spread spectrum system. Some of these 
results are presented in [33]. First the PN spread spectrum communication model is 
presented, followed by the performance of the system with respect to the filtering of stable 
and swept tones. 

4.5.1 PN Spread Spectrum Communication System Model 

The communications model used in this section and shown in Fig. 14 is the 
low—pass equivalent form of an undistorted BPSK signal with additive white Gaussian 
noise. The PN sequence is modelled as rectangular chips of amplitude ± 1 and a chip rate 
R c = 1 Hz. The information bits dk are likewise rectangular, of amplitude ± 1, and of 
duration Tb = L/R c , where L is the number of chips per bit. The signal is corrupted by 
white Gaussian noise, n n , of variance <r nse 2 = N c ,/2, and interference i n of power p(0). The 
effect of the integrate—and-dump circuit on the interference as shown in Fig. 2(b) has not 
been included in the model. 

Two adaptive filter algorithms are considered, i.e., the block least squares and 
recursive least squares algorithms. The block algorithm calculates a set of tap coefficients 
for a transversal filter using the software in [2]. (An alternative would have been to use the 
software in [48],[49]). The recursive algorithm listed in Table 3 calculates a set of forward 
and backward reflection coefficients for a lattice least squares filter, as well as the forward 
and backward prediction errors. 

The filtered data are applied to the PN correlator whose output is summed and 
sampled every L chips. 

For the BLS algorithm, two FIR filter transfer functions will be considered, i.e., 

A(z) and A(z)A'(l/z'), where A'(l/z‘) is the matched filter of A(z). This latter 
configuration causes a delay of p samples, where p is the filter order of A(z). Thus, the PN 
sequence in the correlator must be delayed by p chips for this case. 

As discussed in Section 4.3.3, the RLS algorithm can have associated with it a 
weighting parameter or "forgetting factor" A, where 0 < A < 1 . The effect of this parameter 
on the interference rejection capability of the adaptive filter for stable and swept tones will 
be examined. Finally, the matched filter for the lattice RLS algorithm is obtained by 
feeding the forward prediction error from the i ,h lattice stage into a second lattice having 
the same reflection coefficients. The output is the backward prediction error of the i th 
stage of that filter as shown in Fig. 15. The output is delayed by i samples. (If one was to 
use the FTF least squares algorithm, for example, then the matched filter would be the 
backward predictor.) 
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(C) 

FIGURE 13: SEVERAL ALTERNATIVE DECISION FEEDBACK STRUCTURES 
BASED ON (a) CHIP DECISIONS, (b) IMPROVING BIT DECISIONS 
AND (c) A MATCHED FILTER ARRANGEMENT. 
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FIGURE 14: COMMUNICATIONS MODEL USED IN THE COMPUTER 
SIMULATIONS. 
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FIGURE 15: RECURSIVE LEAST SQUARES LATTICE STRUCTURE FOLLOWED 
BY ITS MATCHED FILTER. 
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4.5.2 


Stable Tone Performance 


The performance of the pre—windowed RLS lattice algorithm is examined with 
respect to a stable tone located at 0.2 Hz. The objective is to evaluate the effect of the 
filter order p and forgetting factor \ on the amount of interference suppression that the 
algorithm can provide. 

The results are presented in terms of: (a) root locations of the whitening filter on 
the z—plane; (b) the average attenuation at the interference frequency; and (c) bit error 
rate performance. To determine the root locations of the whitening filter, the coefficients 
of Am must be calculated from the lattice filter. This can be done using the algorithm 
listed in Table 4 [12],[13] in conjunction with the one listed in Table 3. The algorithm in 
Table 4 is a more generalized version of the LD algorithm, developed for the case of 
pre—windowed signals. All filter orders for A(z) up to and including p are determined. 

The coefficients of each order are roresented as {ao<°)}, {aoO), a^ 1 )}, {aoH a^p), • • •, a p (p)} 
where ao (i ) = 1. Once the coefficient of A(z) are known, the attenuation at the tone 
frequency can be calculated. 

The tone is defined as 

i n — Acos(2T#n + 9 ) (63) 

where A = 1.0, 6i = 0.2 Hz, and 9 is a uniform random variable between 0 and 2 t radians. 
The interference-to-noise ratio is 20 dB. 

Figs. 16 and 17 show the distribution of the roots of A(z) for filter orders ranging 
from p = 2 to p = 8. Corresponding to each filter order is a weighting parameter of either 
1.0 or 0.9. Table 5 lists the average attenuation at the tone frequency achieved for each 
situation. This average was calculated on the basis of 150 sample points, excluding the 
initial transient. The same data were used in each test. 

There are two observations that can be made from Figs. i6 and 17, and the data in 
Table 5: 


(a) for \ = 1.0, the attenuation improves as the filter order increases; it 
increases steadily until p = 6; at p = 8, there is no change. One pair of 
zeroes is dedicated to removing the tone, the others aTe contained within 
the circle and provide little attenuation at the frequencies at which they 
are located as determined by their distance from the unit circle. 

(b) for A = 0.9, the attenuation is less, and does not improve with increasing 
order. The roots wander around more, and in many cases stray towards 
the edge of the unit circle, suggesting that deep notches would be created 
at these extraneous frequencies. 

These observations lead to the conclusion that the RLS algorithm is expected to perform 
better for the case of 1 = 1.0 and a filter order of 6 or 8. This latter point makes sense, 
since a sinusoid in noise is best modelled as an autoregressive moving average ^ARMA) 
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Table 3 


Unnormalized Pre—Windowed Least Squares Lattice Algorithm 

Input parameters: 

p = maximum order of lattice 
x T = data sequence at time T 

X = exponential weighting factor 
Variables: 

R[ T ) R^t-i = sample variances of forward/backward errors 
A iT = sample partial correlation coefficient 

7 ? T = 1 — 7 j T = likelihood variable 

e iT , r iT _[ = forward/backward prediction errors 

K;,t> K*t = forward/backward reflection coefficients 

These computations will be performed once for every time step (T=0,- • • ,T ). 
Initialize: 

® 0 T ~ r 0 ,T“ X T 

Rj x = Rq T = ARJ t .j + X* 

“7 C — 1 

7 -l.T 1 


Repeat for i = 0 to min (p.T) — 1 

A ;<1t = -lA i<l!T .] 4- e iT r iT .j / 7 h.t-i 

7?,t = ~ r i, T / R i:, 

R i\l.T = ^ i♦ l.T / RV.I-, 


e i*l,T “ e i,T R i > «I,T r i,T-l 

R[.1,T = R[,T ~ R 


K 


f 

i*l.T 


= A 


i»i,T 


/R 


f 

i,T 


r i*l,T - r i.T- i Rf*l.T e i.T 

R^l.T = Ri.T-1 ~ ^ MljR i*l,T 
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Table 4 


Computing Predictor Coefficients from the Lattice Parameters 


Initialize: 


b-W = 0 for i = 0, • • • ,p—1 
c 0 <°> = 0 

Outer Loop 

Repeat fori = 0. ■ • • ,p 
a 0 (i) = b 0 {i) =1 i = 0 

a 0 (i) = b 0 !i) = 0 i > 0 

Inner Loop 

Repeat for i = 0. - ■ • ,p—1 

bj (i) = bj (i ) - rjCj (i) /7j-i c 
Cj .l (i) = Cj (i > - rjb/b / Rjb 
aj./b = a/b - Kj,i b bj (i 'b 

bj*, (i) = bp l) - K j(1 f aj (i > 

End of Inner Loop 
End of Outer Loop 


Table 5 


Attenuation Level of Single Tone 



1 

Filter Order 

Attenuation (dB) 


1.0 

2 

47 

0.9 

2 

42 

1.0 

4 

53 

0.9 

4 

42 

1.0 

6 

57 

0.9 

6 

42 

1.0 

8 

57 

0.9 

8 

41 
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FIGURE 16: ROOT LOCATIONS OF WHITENING FILTER WHEN THE 
INTERFERENCE IS A SINGLE TONE; (a) p = 2, A = 1.0; 

(b) p = 2, A = 0.90; (c) p = 4, A = 1.0; 

(d) p = 4, A = 0.90. 














FIGURE 17: ROOT LOCATIONS OF WHITENING FILTER WHEN THE 
INTERFERENCE IS A SINGLE TONE; (a) p = 6, A = 1.0; 

(b) p = 6, A = 0.90; (c) p = 8, A = 1.0; 

(d) p = 8, A = 0.90. 
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process [4],[29], which is equivalent to an AR process of infinite order. The effect of the 
additional roots is to model the background noise (two of the roots are dedicated to the 
tone). For p = 2, the filter must balance its resources to model both signal and noise. 

Thus its notch will not be as deep. As p increases, more resources are allocated to model 
the noise. 

For values of \ less than 1.0, the effective memory of the adaptive filter is reduced, 
so that there is not as much data averaging. This exponential window also leads to 
fluctuations in the filter coefficients, which increase the movement of the roots of the 
rejection filter. 

The above observations for A = 1.0 do not necessarily mean that the excision filter 
will perform better as one goes to higher and higher filter orders. Recall that the objective 
of the excision filter is to whiten the sum of signal, noise and interference. As the filter 
order increases, more signal distortion can be expected and at some point will become a 
significant contributor to noise in Eq. 59(b). It thus becomes a tradeoff between the three 
noise terms in Eq. (59b) as to which one will dominate. 

Bit Error Rate Performance 


This next part deals with the bit error rate performance of the RLS algorithm as 
applied to a stable tone defined in Eq. (63). Both theoretical and simulation results will be 
presented. The model for calculating the theoretical BER is shown in Fig. 18. This model 
is similar to the one in [26], in which the LMS algorithm performance was assessed. This 
type of model is required for the RLS algorithm when 4 < 1.0 which, as mentioned above, 
leads to filter fluctuations and an increase in noise, known as the excess noise. 

The average mean squared value of this excess error is [10] 


= E{Aa'R A a} 

(64a) 

= -nr Tr U + E(p,fl Ko Pl 

(64b) 


where Aa is the vector of filter misadjustments due to the exponential window, R is the 
autocorrelation matrix of the data, Tr refers to the trace operator, £ opt is the optimum 
mean squared error obtained from the Wiener-IIopf solution, I is the identity matrix, and 
P n is a zero mean fluctuation matrix. A simplification [10] of Eq. (64b) is 


E 


excess 


l - i 
1 + X 



1 - 1 
1 + X 



(65) 


where 7 = var(x n 2 )/var 2 (x n ). If x n is a Gaussian process, then 7 = 2 . 

To determine the BER for the system in Fig. 18, the test statistic at the output of 
the bit detector must be determined. This test statistic is 


U — Put f"n “h 1 n. excess) 


n=l 


(66a) 
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FIGURE 18 
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MODEL FOR A LEAST SQUARES EXCISION FILTER FOR THE 
CASE WHEN 1 < 1.0. 
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(66b) 


Li Ij Li Lj 

— Pn^n P n "h Pn- 

n=l n=l i=l n=l 


n.excess 


Letting x n = s n + in + n n , and assuming that the excess error is zero mean and the chips 
are independent samples, then E{U} = L in Eq. (66b). The variance is similar to the 
denominator in Eq. (59b) except for the term due to the excess error. Thus, 

| P P P P P 

var(U) = L< ^aj 2 + /?(0) - 2^\jp(j) 4- ^ ^aja^(j-k) + (r^e 2 + <r nS e 2 ^T aj 2 + £ exce s 

U=1 j=l j =1 k=l j=l 


= <r 2 . (67) 

If the assumption is made that U is a Gaussian random variable, with var(U) = a 2 , then 
the probability of a bit error is defined as 


P e = —-- [ e' ^ dz (68) 

4~TT J 

v/sm; 

where SNR 0 = L 2 / <r 2 . Excess is a function of 7, which can be calculated from Eq. (65). 

For the process under consideration, 

var 2 (x n ) = [1 + (Tn.se 2 +/?(0)] 2 (69a) 

var(x n 2 ) = 2<r nse 2 ((r Ilse 2 + 2) + 4 /j( 0){1 + (T US f 2 } + j p 7 { 0)- (69b) 

The theoretical BER curves as determined from Eq. (68) are shown in Fig. 19. The 
conditions here correspond to a forward prediction filter without its matched filter, and 
filter orders of 2, 4, 6, and 8. The weighting factors were A = 1.0, 0.95, and 0.90. The 
frequency of the interference was offset from the carrier by 0.2 Hz. The results show that 
performance improves as the filter order increases. After close examination of Figs. 19(a), 
(b), and (c), one can see the effect of the excess noise as \ decreases. 

Computer simulations of the communications model depicted in Fig. 14 were also 
carried out. They are based on the following parameters: 

Signal—to—Interference Ratio (SIR): -20 dB 
Processing Gain: 20 

Sampling Rate: 1 Hz 
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FIGURE 19: THEORETICAL BER CURVES BASED ON EQUATION 68 FOR 

SINGLE TONE INTERFERENCE LOCATED AT 0.20 HZ; (a) } = 1.0; 
(b) \ = 0.95; (c) 1 = 0.90. 
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FIGURE 20: BER CURVES BASED ON SIMULATION FOR SINGLE TONE 
INTERFERENCE LOCATED AT 0.20 HZ; (a) A = 1.0; (b) A = 0.95; (c) A = 0.50. 















Slightly more than 25000 random bits of information were generated, yielding over 
500 thousand chips after multiplication by a random PN sequence. To each signal sample 
(chip) was added a sample of white Gaussian noise n n and interference i n - The noise 
variance was determined from the Eb/N 0 of interest, i.e., a nse 7 = L/(2(Eb/N 0 ) where N 0 is 
the single—sided noise power spectral density. The amplitude of the interference was set 
equal to 14.142 to provide an SIR/chip of —20 dB. The BER was determined by comparing 
the bits "transmitted" to those "received". 

The results are shown in Fig. 20. ForA = 1.0 and 0.95 [Figs. 20(a) and (b)] there is 
good agreement between theory and simulation. However, for } = 0.90 (Fig. 20 (c)), the 
simulated results are more pessimistic, indicating that the approximation to Eq. (64b) 
(namely Eq. (65)) may be invalid for such low values of One would therefore have to 
use the more general Eq. (64b) for a more accurate estimate of the excess error. 

The next section considers a swept tone interferer. The performance of the BLS and 
RLS algorithms are presented in terms of their tracking ability and the bit error rate at the 
output of the bit detector. The results were obtained from computer simulation. 

4.5.3 Swept Tone Performance 

The interference is modelled as a frequency modulated waveform of the type 

t 

i(t) = Acos[2T<5ft + d f m(r)dr + 0] (70) 

0 

where Sf is an offset frequency from the spread spectrum signal carrier, m(r) is the 
modulating signal, d is the frequency deviation constant of the modulator in 
radians/second/volt [85], and 0 is an initial random phase. 

To demonstrate the tracking behaviour of the BLS and RLS adaptive filters, a 
swept tone with forward and backward sweep rates of 0.001 Hz/sec was modelled using Eq. 
(70). A section of m(r) is shown in Fig. 21(a), sweeping between 0.1 and 0.3 Hz, and 6f is 
0.2 Hz. The amplitude of m(r) was set to unity and d was chosen to be 0.628 
radians/second/volt. The peak frequency deviation Af for this case was 0.1 Hz. The 
interference—to—noise ratio was set at 20 dB. The spread spectrum signal was not included 
in this set of tests. 

The BLS method was examined first. For each block of data, the set of 
autoregressive parameters {a,, i = 1,2, * • -,p, £r p 2 } were obtained along with the roots of A(z), 
which represent the locations of the notches created by the filter. 

Figures 21(b) to 21(d) show three examples of the tracking capability of a fourth 
order BLS filter A(z). The number of samples per block is, respectively, 100, 50, and 10. 
Only the positive frequency is plotted. The horizontal lines correspond to the frequencies 
of two of the roots of A(z) in the upper half of the complex z—plane, the other two roots 
being in the lower half plane. The location of the horizontal lines indicates the position of 
the notches for the particular block of data processed to produce the adaptive filter 
coefficients. 
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FIGURE 21: TRACKING RESULTS FOR THE BLS ALGORITHM FOR p = 4 AND 
VARIOUS BLOCK SIZES N, (a) MODULATING SIGNAL OF FM 
INTERFERER; (b) N = 100, (c) N = 50, (d) N = 10. 
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FIGURE 22: TRACKING PERFORMANCE OF THE RLS ALGORITHM; (a) p = 2, 

A = 0.90; (b) p = 4, A = 0.90; (c) p = 2, A = 0 95. 
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(a) 


(b) 


FIGURE 23: TRACKING PERFORMANCE OF THE RLS ALGORITHM; (a) p = 4, 
A = 0.99; (b) p = 4, A = 1.0. 
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For the case of 100 samples per block, Fig. 21(b) shows several instances when the 
notch filter is outside the range of the tone. One would therefore expect poor tone 
attenuation during these times. Reducing the number of samples per block to 50 (Fig. 
21(c)) improves the tracking capability somewhat, but at times only two out of the four 
roots follow the tone. Fig. 21(d) shows good tracking of the tone when only 10 samples per 
block are used. In this case, one pair of roots track the tone, whereas the other two are 
located at 0.0 and 0.5 Hz most of the time, i.e., on the real axis of the complex z-plane. 

For the RLS lattice algorithm, the samples were processed recursively using the 
prewindowed lattice algorithm listed in Table 3. The filter parameters of A(z) for this 
lattice were obtained from the algorithm listed in Table 4. 

Figures 22(a) and (b) show the tracking capability of the RLS lattice filter for the 
case of A = 0.90 ana filter orders of 2 and 4, respectively. As in the BLS case, for a fourth 
order filter, two of the roots in Fig. 22(b) track the tone, whereas the other pair wander 
around, similar to the case for the stable tone. 

These results lead to some preliminary conclusions concerning the BLS and RLS 
filters. First, in the case of the BLS filter and a fairly fast sweeping tone covering a large 
percentage of the spread bandwidth (in this example, the swept tone covered a bandwidth 
equal to 20% of the chip rate), as few samples as possible per block should be used to 
achieve reasonable tracking and filtering. Second, for both the BLS and RLS algorithms, a 
second order filter appears to be more appropriate, since the additional roots from higher 
order filters have a tendency either to locate themselves on the real axis or elsewhere. 

These extra roots will introduce additional spectral distortion in the spread spectrum 
signal. The degree of distortion would depend on the depth and quantity of the notches. 

As demonstrated for the stable tone interferes for decreasing values of A, the roots were 
less confined to one spot in the interior of the unit circle and had a tendency to wander 
quite significantly. The notch depths would be enhanced if the higher order filter was 
cascaded by its matched filter A'(l/z’). For larger values of A, the tracking capability of 
the filter would degrade, i.e., a larger lag would be introduced [Fig. 22(c)]. Thus, the 
amount of interference attenuation achieved would be less, thereby increasing the BER. As 
A approaches 1.0 all the roots contribute to the filtering; they also exhibit much less 
mobility [Fig. 23(a)]. In the limit, when A = 1.0 the roots tend to become stationary [Fig. 
23(b)], In other words, the RLS algorithm approaches the Wiener—Hopf solution. 

Bit Error Rate Performance 


The simulation conditions were the same as those used for the stable tone. Figs. 
24(a) and (b) show the BER’s for the BLS method for second order (p = 2) and fourth 
order (p = 4) filters A(z), respectively; the matched filter results are also included here. In 
these and subsequent figures, the BER when the excision filter was not present is not 
shown; however, for this case it was approximately 0.4 for all of the Eb/N 0 ’s considered. 

In general, the BER improves as the block size decreases, confirming what was 
observed previously when the BLS algorithm’s tracking capability was explored (Fig. 21). 
First consider the results for N = 100. When the matched filter was not used, there is a 
slight improvement in performance in going from p = 2 to p = 4, suggesting that the 
increased bandwidth provided by the higher order filter suppressed more of the interference 


- 57 - 



Eh/No (db) 


V^ovdB) 


FIGURE 24: BER RESULTS OF BLS AND RLS ALGORITHMS WHEN THE 
INTERFERENCE IS A SWEPT TONE. THE DASHED LINE REFERS TO 

THE CASES IN WHICH THE MATCHED FILTER WAS 
INCLUDED; SIR/CHIP = -20 dB. 

(a) BLS ALGORITHM, p = 2, "A”, "X", AND "0" REFER TO 
N = 100, 50, AND 10 SAMPLES/BLOCK, RESPECTIVELY; 

(b) BLS ALGORITHM, p = 4, N = 100, 50, AND 10 
SAMPLES/BLOCK; 

(c) RLS ALGORITHM, A = 0.90. p = 2 ("O") AND p = 4 ("X"V 

(d) RLS ALGORITHM, p = 2, A. = 0.90 ("O”), A = 0.95 ("X"). 






(note that from Fig. 21 all the roots took part in the filtering process). When the matched 
filter was included, the performance improved significantly, particularly when the filter 
order was increased to p = 4, indicating a much larger reduction in the interference for this 
case. 


For N = 50 and no matched filter, the BER results reveal very little difference 
between the second order and fourth order filters. Referring to Fig. 21(b), we see that for p 
= 4 there is a tendency for two of the roots to track the tone and the other two to 
sometimes wander around. Furthermore, because of the almost imperceptible difference 
between the results for p = 2 and p = 4, one comes to the conclusion that the wandering 
roots do not introduce much signal distortion, that the noise in Eq. (59b) is dominated by 
the residual interference. However, the matched filter results indicate slightly inferior 
performance for the fourth order filter, suggesting that perhaps signal distortion is starting 
to manifest itself. For N = 10, this effect is more evident. 

Figures 24(c) and (d) illustrate the BER performance for the RLS lattice filter. The 
performance of both the second and fourth order filters for A = 0.95 are shown in Fig. 

24(c). Consider first the case without the matched filter. The results exhibit a 
degradation in performance for the fourth order filter, indicating an increase in the noise at 
the output of the adaptive filter. This noise is likely dominated by the residual 
interference. When the matched filter was present, the performance between the two cases 
was very similar, suggesting two things: (a) they both removed about the same level of 
interference; and, (b) a filter order of p = 4 did not appear to increase the other noise 
terms (signal distortion, Gaussian noise, and the excess noise due to filter coefficient 
fluctuation). More tests using higher order filters would have to be conducted to examine 
this problem in more depth. 

Similar comments hold for the results in Fig. 24(d) for the case of a second order 
filter, and A = 0.90 and 0.95. Recall from Fig. 22 that increasing A resulted in an increase 
in the lag between the predictor and the interferer’s modulating signal. This would result 
in an increase in the residual interference and, therefore, an increase in BER as shown 
when the matched filter was not present. With the matched filter, the effect of the lag is 
not as noticeable. 


5.0 CONCLUDING REMARKS 

This report has presented a review of the research which has been conducted in the 
area of jammer suppression in DS/SS systems. Most of the emphasis has been on digital 
adaptive schemes employing linear and non-linear (i.e., decision-feedback configurations) 
filtering. Furthermore, the majority of effort had been on the filtering of single or multiple 
tones, as well as narrowband coloured noise. Very little work had been done on 
interference of the FM-type. The latter part of this report shed some light on the 
problems associated with filtering a swept tone interferer. 

The results suggest that for a single swept tone, a second order filter will achieve 
better performance than higher order filters, since the latter may introduce additional 
signal distortion. This was particularly evident with the BLS algorithm and small block 
sizes. For the RLS algorithm, it was shown too that a second order filter is adequate 
enough for the suppression of this type of interference; signal distortion did not appear to 
present much of a problem for a fourth order filter. One would expect a degradation in 
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performance for higher filter orders and X < 1. To be more conclusive, additional work 
would have to be conducted in this area. Also indicated was the sensitivity of BER to 
changes in A, especially when the matched filter was not used. Finally, the results 
produced by the BLS and RLS algorithms contrast those of the stable tone in which 
increasing the filter order improved the performance. 

The tracking examples involving the RLS algorithm exhibited a lag between the 
modulating signal and the e:timate. This lag is a function of A and suggests that an 
artificial delay could be imposed on the data to compensate for it. One would therefore 
expect an improvement in performance if this were done. 

From the results presented for the swept tone, the block and recursive least squares 
second order filters can be viewed as tunable notch filters or, conversely, as frequency 
demodulators if one calculates the roots of the excision filter over time. This concept opens 
the door to other possible filtering schemes. 

With this in mind, the author has developed a new digital adaptive algorithm which 
is based on an extended Kalman filter [68], Kelly and Gupta [30] and Polk and Gupta [63], 
using Kalman filtering theory, developed a digital PLL which has performance 
characteristics similar to that of an analog PLL. This model forms the basis of a new 
excision algorithm to be described in a subsequent report. Their model has been modified, 
and extended to render it adaptive. In particular, the performance of the new algorithm 
will be presented for the case of FM interferers such as the type considered in this report. 
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6.0 GLOSSARY 


AR - Autoregressive 

BER — Bit Error Rate 

BLS - Block Least Squares 

BPSK — Binary Phase Shift Keyed 

BW — Bandwidth 

CCD — Charge Coupled Device 

CTD - Charge Transfer Device 

DPT — Discrete Fourier Transform 

DS/SS — Direct Sequence Spread Spectrum 

E{U} — Statistical Expectation of the Random Variable U 

Eb/No — Energy per Bit—to—Single-Sided Noise Power Spectral Density Ratio 

ESM — Electronic Support Measures 

FAEST — Fast a posteriori Error Sequential Technique 

FFT — Fast Fourier Transform 

FIR — Finite Impulse Response 

FM — Frequency Modulation 

FTF — Fast Transversal Filters 

IF — Intermediate Frequency 

I/S — Interference—to—Signal Ratio 

LD — Levinson—Durbin 

LMS - Least Mean Square 

LPF — Low Pass Filter 

LS - Least Squares 

Ml —Mixer 

M2 — Mixer 

ML — Maximum Likelihood 

PLL — Pha^e-Locked Loop 

PN — Pseudo Noise 

QPSK — Quadrature Phase Shift Keyed 

RLS — Recursive Least Squares 

SAW — Surface Acoustic Wave 
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